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Abstract. The nonequilibrium stationary state of an exclusive genetic switch is 
considered. The model comprises two competing species and a single binding site 
which, when bound to by a protein of one species, causes the other species to be 
repressed. The model may be thought of as a minimal model of the power struggle 
between two competing parties. Exact solutions are given for the limits of vanishing 
binding/unbinding rates and infinite binding/unbinding rates. A mean field theory is 
introduced which is exact in the limit of vanishing binding/unbinding rates. The mean 
field theory and numerical simulations reveal that generically bistability occurs and 
the system is in a symmetry broken state. An exact perturbative solution which in 
principle allows the nonequilibrium stationary state to be computed is also developed 
and computed to first and second order. 
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1. Introduction 

Genetic networks are interacting, many-component systems of genes, RNA and proteins, 
that control the functions of hving cells [1-3]. They exhibit complex, nonlinear 
behaviour due to the feedback between the expression of different genes. At a simplistic 
level of description the gene sequence, through transcription by RNA and translation 
by mRNA, leads to the production of proteins. These proteins may themselves act 
as transcription factors which may switch on or off the activity of other genes [4]. A 
particularly simple example of a genetic network is a toggle switch formed from pairs 
of genes that mutually repress each other's expressions [5]. Such switches serve as 
microscopic models for bistable and oscillatory states [6-9] and importantly may be 
synthesized [10]. 

Genetic switches may exhibit bistability where there are two possible dynamically 
stable long-lived states for the switch. There has been considerable interest in how 
bistability may be maintained and how it is effected by stochastic fluctuations due to 
small numbers of proteins and intrinsic noise [11-13]. In particular the switching time 
between the two states has been measured and numerical techniques have been devised 
to study the switching time within theoretical models [7, 14-16]. 

In this paper we consider the simplest toggle switch introduced by Warren and P.R. 
ten Wolde [6] which we will refer to as the Exclusive Switch (this model is referred to 
in [7] as the exclusive switch without cooperative binding). The model comprises two 
genes labelled 1,2 each leading to the production of proteins Xi and X2. These proteins 
also degrade stochastically. There is a single mutual binding site to which either an Xi 
or X2 may bind and when bound repress the production of the other protein. Thus when 
an Xi is bound, the Xi population fluctuates around some steady state value determined 
by the balance of production and degradation while the X2 population degrades towards 
zero. 

This switch may be considered more generally as a minimal model for the power 
struggle between two competing parties. This is illustrated by the following charicature 
of "mob dynamics" . There are two competing parties or gangs of individuals and room 
for only one individual to wield absolute power. When an individual of one party is in 
power his own party membership may grow but the other party membership dwindles. 
Random influences imply that the control of power is occasionally lost and power may 
be seized by any individual. Thus in a temporary power vacuum the membership of the 
minority party will increase and there is a small chance that a member of this party will 
seize power and that the minority will eventually become the majority. 

In the context of statistical physics the Exclusive Switch is an example of a 
nonequilibrium system. This is because the microscopic stochastic dynamics do not 
obey detailed balance. For example when an Xi protein is bound the degradation of 
an X2 protein is irreversible. The structure of nonequlibrium stationary states has been 
of considerable interest and is generally characterised by the existence of probability 
currents in the stationary state (which do not exist in equilibrium stationary states 
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due to the presence of detailed balance) [17]. For example it has been shown that 
spontaneous symmetry breaking may occur in non-equilibrium systems under conditions 
where it is precluded from their equilibrium counterparts e.g. in one spatial dimension 
[18,19]. 

The bistability exhibited in the Exclusive Switch may be thought of as symmetry 
breaking where although the microscopic dynamics is symmetric between the two 
proteins, the stationary state comprises two possible long-lived dynamical states in 
which the symmetry is broken and one protein dominates. In an equilibrium system 
the switching time between the two symmetry- broken states may be estimated by the 
Arrhenius law r ^ exp f3AF where the free energy barrier AF is extensive in the 
system size. For a nonequilibrium system on the other hand the free energy or indeed 
the stationary state is not known a priori and one is required to construct the stationary 
state on a model by model basis. For genetic switches there has been recent interest in 
developing analytical approaches to describe the stationary states [20,21]. 

In the present work we study analytically the nonequilibrium stationary state of the 
Exclusive Switch. Previous analytical studies have concentrated on systems with only 
one gene [20,22,23]. Our aim is to understand whether symmetry breaking occurs and, 
if so, the nature of the symmetry broken state. To this end we develop two analytical 
approaches. First we construct a mean field theory. Then we develop a perturbative 
approach that in principle allows the nonequilibrium stationary state to be computed 
exactly and we present analytical results to first order. A complementary approach 
to this system has been developed in [24], where an approximation scheme based on 
eflFective interactions is used. 

The paper is organised as follows. In section 2 we define the Exclusive Switch model 
and write down the system of master equations that describe the system. We also present 
some numerical simulations which illustrate the nature of the symmetry breaking and 
consider exactly solvable limits. In section 3 we present a mean field theory and compare 
to stochastic simulations. In section 4 we develop a perturbative approach, compute 
the results to first order and compare with simulation results. Conclusions are drawn 
in section 5. 

2. Model Definition 

The state of the system is defined by: the number of free proteins of type 1, A^i ; the 
number of free proteins of type 2, A^2, and the state of the switch, S', which takes value 
if no protein is bound, value 1 if a protein of type 1 is bound, and value 2 if a protein 
of type 2 is bound. 

The stochastic dynamical processes are as follows: a protein degrades (leaves the 
system) with rate rf; when the switch state is proteins of both type 1 and 2 are 
produced with rate g; when the switch state is 1 proteins of type 1 produced with rate 
g and when the switch state is 2 proteins of type 2 produced with rate g] if the switch 
state is a protein binds with rate b and when binding occurs the switch state changes 
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to the type of bound protein and the number of free proteins is reduced by 1; a bound 
protein unbinds with rate u and when unbinding occurs the switch state changes to 
and the number of free proteins is increased by 1. 

We shaU consider the joint probabihties Ps{Ni^ N2) of the protein numbers A^i, A^2, 
switch state S. 

2.1. Master equation 

FoUowing the stochastic dynamical processes described above the system of master 
equations that defines the model can be written as follows 

= ^[Po(iVi - 1, iV2) + Po(iVi, iV2 - 1) - 2Po(iVi, iVs)] 
+ d[(iVi + l)Po(iVi + 1, N2) + {N2 + l)Po(iVi, N2 + 1) 

- (iVi + iV2)Po(iVi, N2)] - b{Ni + N2)Po{m, N2) 
+ u[Pi{N, - 1, N2) + P2{Ni, N2 - 1)] (1) 

= g[P,{Ni-l,N2)-Pi{N,,N2)] 
+ d[{N^ + l)Pi(A^i + 1, A^2) + (A^2 + N2 + 1) 

-{Ni + N2)Pi{N,,N2)] 
+ b{N, + l)Po{Ni + 1, N2) - uPiim, N2) (2) 

= g[P2{Ni,N2-l)-P2{N,,N2)] 

+ d[{N^ + l)P2(A^i + 1, N2) + {N2 + l)P2(iVi, N2 + 1) 

-{Ni + N2)P2{m,N2)] 

+ b{N2 + l)Po(iVi, N2 + I)- uP2{m, N2) (3) 

where g is the generation rate of a protein (when the generation is not suppressed by 
the switch state), d is the degeneration rate of a single protein, b is the binding rate of a 
single protein and u is the unbinding rate of the bound protein. The whole problem is 
clearly symmetric with respect to the variables 1 and 2. Also note that, the degeneration 
term is the same for the three equations, since it does not depend on the value of 5'. 

2.2. Nature of symmetry breaking 

In order to illustrate the qualitative behaviour, we first present stochastic simulations 
of the Exclusive Switch. These simulations are performed with a Gillespie algorithm 
[25,26], in which the reactions described in the model are given a certain probability, 
depending on the state of the system and the value of the parameters. These 
probabilities are used to determine stochastically which reaction is going to happen 
next, and when it will happen. The time that the system spends in a given state 
A^i, iV2, S' is normalized to obtain the probability distributions. The reference values for 
the simulations are the typical values for bacteria such as Escherichia coli [1] 

^ = 0.05, rf = 0.005, 6 = 0.1, « = 0.005 (4) 



dPo 
dt 
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dP2 
dt 



iNi,N2) 
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In the system there is a clear symmetry between the two proteins species, since 
they undergo the same microscopic reactions. However, at any given time the system 
is typicaUy dominated by one of the proteins; thus the symmetry is broken. The reason 
for this is that, once a protein, e.g. of type Xi, binds to the promoter site, proteins 
X2 start disappearing, while the number of Xi fluctuates around a steady value. That 
means that when the bound protein unbinds (as it will do eventually due to stochasticity 
of the system), it is much more probable for proteins Xi to bind to the promoter site 
again, since there are more of them. At the same time it is more diflicult for proteins 
X2 to bind to the promoter site. However the X2 will not disappear permanently from 
the system (there is no absorbing state), and will be produced again the moment the 
bound protein Xi unbinds. Thus, with a small probability, proteins X2 will be able to 
bind again to the promoter site, and become the dominant species, as proteins Xi start 
to degenerate. This means that there are two symmetry-broken states, in which one 
species is much more abundant than the other. 

With regard to the probability distribution P(A^i, A^2), this bistability is translated 
into two peaks, concentrated around the axes, i.e, where one of the protein numbers is 
almost zero. This is illustrated in a contour plot in flgure 1 for u = 0.05. However, 
this bistability depends strongly on the value of the parameters: the bigger the value of 
the more irrelevant is the switch state for the dynamics of the protein, and the less 
important is the bistability we have described. For example in flgure 1, when g^b and d 
are kept constant and u is increased, the peaks move together and eventually merge at 
some value oi u ^ 0.15. Thus there is an apparent transition from a distribution with 
two symmetry-related peaks to a distribution with one symmetric peak. We wish to 
study the nature of this transition i.e. is there an underlying phase transition at a flnite 
value of u where the system changes from symmetry- broken behaviour to symmetric 
behaviour, or is the transition simply due to two peaks coming closer together and no 
longer being resolved? The latter would correspond to a 'geometrical transition' in the 
form of P{Ni^ N2) but would not correspond to any underlying phase transition. 

Although the whole probability distribution is always symmetric, distributions Pi 
and P2 will be a priori asymmetric, since they describe the probability of the number 
of proteins when a protein 1 or 2 are bound, which are not symmetric situations. Let 
us now deflne ta and as the probabilty masses of Pi on either sides of the diagonal 

1 (5) 
rB= ^1(^1^ ^2) + - ^i(^i'^2). 

We can now study how Ta^Tb change with and whether there is a clear transition 
between the situation in which they are different, and the one in which they are equal 
to each other (if there is any). Figure 2 shows that these two quantities approach each 
other in a continuous way, and they are equal only when 2/ ^ oc, that is, when the 
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u=0.05 u=0.1 




Figure 1. Contour plot of probability distribution P(A^i, A^2) obtained from stochastic 
simulations. A Transition from a two peak regime to a one peak regime occurs as the 
unbinding parameter u increases, g^b and d are kept equal to the E. coli values (4) 

only possible state of the switch is S' = and it has no longer any effect on the protein 
dynamics. Therefore, the probability distributions Pi, P2 and hence ta^Tb tend to zero, 
all the probability being concentrated in the distribution Pq (^15^2)5 which is always 
symmetric. 

Despite the fact that a change in the shape of the distribution is observed (figure 
3), there is no evidence of the typical singularities that appear in a phase transition. 
Instead Pi appears to deform continuously into a symmetric distribution when u ^ 00. 
We conclude that (at least for these parameter values) Pi(A^i, A^2) remains asymmetric 
for 2i < oc, and as a consequence so does P2, and there is no transition to a symmetric 
state in this marginal distribution. As figure 3 shows. Pi has only one peak for different 
values of and even if it becomes smaller as u increases, it does not become symmetric 
at any point. 

We deduce that, even though P(A^i,A^2) appears to become symmetric at some 
finite value of u (see figure 1), there is no true phase transition between symmetric and 
asymmetric regimes, since Pi and P2 always remain asymmetric. Thus, the bistability 
of the switch is always present, with the asymmetric distributions Pi,P2 decreasing as 
the switch state becomes less important, that is, as 2/ increases. 
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Figure 2. Evolution of the two contributions to Pi(A^i,A^2) defined in (5) — for 
Ni > N2 and for Ni < N2 — as u increases. 

r — r 




2.3. Solution in limit b ^ 

We now consider the system in the hmit b ^ with the binding constant 

k = - (6) 
u 

held fixed. In this limit the system will equilibrate between binding/unbinding events. 
When the switch is in state 1 the number of type 2 proteins decays to zero and the 
number of type 1 is given by Poissonian statistics for the generation/degeneration 
processes: 

PM,N2) = [^Y' (7) 

where ri is the probability of the switch being in state 1 (ri = ^^^q Pi(A^i, A^2))- 
Similarly, in state 2 the number of type 2 is given by Poissonian statistics 

P,{N^, N,) = r,^ (^Y' e'^/^.o (8) 



and in state both A^i and N2 are given by Poissonian statistics 

1 ^JA"'*"' 



In order to fix the probabilities ro,ri,r2 we consider the master equation for ri which 
can be obtained from summing (2) over the variables Ni^ N2: 

r'l = -uri + bl^ (10) 
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u=0.005 u=0.05 




3 10 13 20 3 10 13 20 



Ni Ml 

Figure 3. Contour plots for the probability distribution Pi(A^i,A^2), for different 
values of u. As u increases, the peak of the distribution moves towards the diagonal 
Ni = A^2, but the distribution is not completely symmetric as long as u is finite. Also, 
the probability mass of the distribution decays as u increases, and is transferred to the 
symmetrical distribution Po{Ni^ N2). 



where the quantity is the mean value of A^i given the system is in switch state 
S = 0^ multiphed by the probabihty of being in switch state S = 0^ i.e, = 
N2=o ^1^1(^15 ^2)- In the stationary state we have = r^g/d and (10) becomes 

-UTi + . (11) 
A similar equation holds for r2 and the normalisation of probability yields 

"° - 1 + 2kg /d l + 2kg/d- ^^^> 

Thus we see that the system has three long-lived states: when S = 1 the system is 
dominated by type 1 proteins; when S = 2 the system is dominated by type 2 proteins, 
and when S = the system is in a symmetric state. 

When k ^ 00 the symmetric state has zero weight {tq 0) and the system is either 
in the S = 1 state or the S = 2 state. The switching time between the two states is 
expected to diverge in this limit therefore the system exhibits symmetry breaking. Note 
that the transition from the S = 1 state to the S = 2 state is through the symmetric 
S' = state. 
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2.4- Solution in limit b ^ oo 

In this limit the unbinding and binding events happen on a faster time scale than the 
growth and degradation of proteins. Therefore the switch state becomes decoupled from 
the numbers of protein and the probabilities obey 

Ps{NuN2) = P{m,N2)rs{Ni,N2) (13) 

where P(A^i, A^2) is the probability of A^i,A^2 proteins being in the system regardless of 
switch state. That is, for given number of proteins A^i and N2 (which include any bound 
protein) the switch probabilities equilibrate and obey 

un{N,, N2) = bNMNu N2) ur^iNr, N^) = bN^r^iN,, N^) (14) 

which may be solved to obtain 

ro = [1 + k{Ni + N2)r^ ri = A:A^iro = kN2ro (15) 

where the binding constant k is given by (6). 

A master equation for the evolution of P{Ni, N2) on the slower timescale on which 
generation and degeneration events occur may then be written down: 
dP 

^(iVi, iV2) = ^[1 - r2(iVi - 1, N2)]P{Ni - 1, N2) - d[Ni - N2)]P{N,, N2) (16) 

- g[l - r2{N,, N2)P{Nu N2)] + d[N, + 1 - n{N, + 1, N2)]P{Ni + 1, N2) 
+ g[l - nim, N2 - l)]P(iVi, N2-I)- d[N2 - r2{m, N2)]P{Nu N2) 

- g[l - ri{Ni,N2)]P{m, N2) + d[N2 + 1 - r2{m, N2 + l)]P{m, N2 + 1) . 

The terms with coefficient g in (16) represent generation of A^^ when the switch is not in 
switch state i. The terms with coefficient d in (16) represent reduction of A^^ with rate 
dNi when the switch is not state i and reduction with rate d{Ni — 1) when the switch 
is in state i (this is because the bound protein does not degrade). The stationary state 
of (16) obeys detailed balance with respect to generation and degradation of A^i and N2 
individually, thus 

= P(N^-l,N2)^- [1 + ^(^1-1)] 1 [l + k{m + N2)] 



d[l + k{Ni + N2- 1)] A^i [1 + A:(A^i + A^2 - 1)] 
with a similar equation holding for detailed balance in A^2- Then these equations may 
be iterated to obtain 

P(N N)- m^^(i±W+M 1 W ^1±^!^ pro AT 1 

P{Nr,N2)- [-) ^ ^ 11 l + fc(iV2 + nO^^°'^^^ 

^ ^ 711=0 ^ ^ 

g^N,+N, (1 + k{Nr + N2)) Un^^ii^ + knr) nS(l + fcn^) 

~d) — Nmi nS"^-^(i + fc-) ^ ' ^ ■ ^ ^ 

The constant P(0, 0) will be determined by normalisation of the sum of probabilities to 
one. 
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The expression (18) is a single-peaked distribution which is symmetric in A^i and N2. 
To see this one can identify the stationary points of P{Ni, N2) through the conditions 
P{Ni + 1,N2) = P{Ni, N2) and P{Ni, N2 + 1) = P{Ni, N2) which yield 

g 1 {1 + kNi) {1 + k{Ni + N2 + 1)) 



diNi + l){l + k{Ni + N2)) {l + kiNi+N2)) 
g 1 {l + kN2) {1 + k{Ni + N2 + 1)) 



1 

1 . 



(19) 



d{N2 + l)il + kiN^ + N2)) il + k{Ni + N2)) 
The solution of these equations is Ni = N2 = N/2 where, when N is large, 



N-2- = . (20) 
a 



Thus, in the u^b ^ oo hmit the system has reached a symmetric state with the 
probabihty distribution peaked at A^i = A^2 = N/2. 

3. Mean Field Theory 

In this section we develop a mean field theory in which some correlations in the numbers 
of proteins are ignored. 

3.1. Exact Moment Equations 

We start from the exact master equations (1-3) for the evolution of probabilities 
Ps{Ni^ N2). The zeroth moments of A^^ are the probabilities rs i.e. 

00 

rs= Yl Ps{NuN2) for 5 = 0,1,2. (21) 

Ni=0,N2=0 

We now define the first moments of Nf of A^^ as follows 

oo 

iVf= NiPs{NuN2) for ^ = 1,2 5 = 0,1,2 (22) 

Ni=0,N2=0 

and the second moments 

oo 

{NiNjf = Y NiNjPs{N^,N2) for i,i = l,2 5 = 0,1,2.(23) 

Summing (1) gives 

^ = _ 5(iVf + iVf ) + u[n + rs] (24) 

Note that the physical meaning of e.g. Nf is the probability of being in switch state S 
(rs) multiplied by the mean number of type i, given that the switch is in state S. 



gro - dN^ - b[{NiNif + (A^iA^s)"] + «[(A^i)^ + {Nif + n] . (25) 
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Similarly, summing (2) gives 

—1 = +m^-uri (26) 
dNl 



dt 



gn - dNl + b{Ni{Ni - 1))° - «(iVi)i (27) 



^ = _ rfATi + biN.N^f - uiN^y . (28) 
We now invoke symmetry between switch state 1 and 2 : 

ri=r2 = (l-ro)/2 =nI = nI ]vf=A^. (29) 
Then the exact steady-state versions of equations (24-28) read 

ri = -M ro = 1 - 2ri (30) 
u 



b[{NiN^y + (iViiV2)o] = gro - dTVO + u[{Niy + {N,y + n] (31) 

b{NiNiy = -gn + {d + u)Nl + bjNif (32) 

b{NiN2f = {d + u)nI . (33) 

Note that if we sum (31-33) we obtain the exact relation 

d{Nf +Nl + Nl)=g{l- n) = g{l - ra) (34) 

which simply gives the overall birth/death balance for A^i. Also note that (32,33) give 
exact relations between the second moments and first moments. However to actually 
evaluate these moments one would have to consider equations for higher moments, 
leading to a hierarchy of equations. 

3.2. Mean Field Approximation 

We now make a mean-field approximation that expresses second moments in terms of 
first moments: 

Tvo Tvo 

(iViiVi)o = — ^ + (35) 

^0 



2 



{NrN^f = = . (36) 

ro ro 

Note that a symmetry condition from (29) has been explicitly used in the last equation 
of (36). The first relation (35) comes from the assumption that A^i has a Poisson 
distribution when the switch is in the state. This means that the second moment is 
equal to the square of the mean plus the mean itself. However, there is an important 
factor To, which comes from the fact that (A^iA^i)^/ro is the mean square value of A^i 
given that the switch is in the state and N^/ro is the mean value of A^i given that the 
switch is in the state. The Poisson approximation is in fact exact in the limit where 
u^b tend to zero (see section 2.3). 

The second relation (36) is a simple factorization scheme which ignores correlations 
between the values of A^i and N2 when the switch state is 0. 
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Using this approximation scheme (32) becomes 



and (33) becomes 



b {N^Y 



gri 



bro 



d + u To d + u d + u 




(37) 




(38) 



Using expressions (37) and (38) in (31), combined with the previous approximations 
(35) and (36) yields the foUowing quadratic equation for N^/tq: 

2 

hg 



2bd 
d + u 




d 



d 



u 




9 = 0. 



One must take the positive root of this quadratic which yields 

]v° 1 r 

—L = — fyg - d{d + u) + {{d{d + u) - hgf + %hdg{d + u)) 
Then using (30) one obtains 



1/2 



1 + 



u ro 




2b m 
1 + 

u ro 



(39) 



(40) 



(41) 



(42) 



One may check the hmits of section 2 from the quadratic (39). In the hmit b^u 
with k = b/u one obtains N^/vq = g/d and Tq = [l + ^] in agreement with Section 
2.3 where it is shown that A^i foUows a Poisson distribution with mean g/d when the 
switch is in state S' = 1 or S' = 0. 

In the hmit b^u oo with k = b/u fixed the quadratic (39) reduces to 

2 




2kd+—(d-kg) - g = . 



(43) 



This quadratic for — is the same as the quadratic (20) for the value of = 2Ni that 



maximises P(A^i, A^2) in the exact solution of section 2.4. 



3.3. Comparison to simulation results 

The mean field theory we have developed can be compared with the simulations 
by studying the zeroth and first order moments of the probabilty distributions, i.e., 
ro, A^i , A^i , A^f . (The probabilities ri and r2 may automatically be obtained from tq, 
since ri = r2 = ^^.) 

DiflFerent values of the mentioned quantities are given in table 1, where the 
parameters of the models have diflFerent values. The reference for this table is the set of 
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E. coli 


MFT 


u=0.05 


MFT 


i 


(4.9186 ±0.0007)- 10-3 
(2.489 ±0.005)- 10-2 
4.942 ± 0.008 
(5.908 ±0.007)- 10-2 


4.92705- 10-3 
2.48768- 10-2 
3.74372 
1.25604 


(4.5263 ±0.0008)- 10-2 
0.23862 ± 0.00014 
4.113 ±0.003 
0.8734 ± 0.0005 


4.54635- 10-2 
0.238634 
2.71128 
2.2774 




u = 50,b= 1000 


MFT 


^i = 5-10-^6= 10-*5 


MFT 


m 


(4.941 ±0.004)- 10-3 
(2.55 ±0.06)- 10-2 
9.89 ±0.11 
0.227 ±0.013 


4.95073- 10-3 
2.48762- 10-2 
2.50019 
2.49969 


(2.48 ±0.03)- 10-3 
(2.477 ±0.024)- 10-2 
4.92 ± 0.05 
(5.012 ±0.007)- 10-5 


2.49872- 10-3 
2.49375- 10-2 
4.98751 
4.97754- 10-5 



Table 1. Results from different quantities from simulations and mean field theory 
approach. The values of ro and predicted by the mean field theory are always in 
good agreement with the simulations. and predictions are quite far from the 
simulation results, except on the limit u^b 0, where our approximation is exact. 



E. coli values, and only the parameters that are changed are written. For the simulations 
performed, ro is always in good agreement with the mean field theory approximation, 
and so are ri and r2. Also, is in quite good agreement, too. 

A^^^ and are different in the mean field theory, which is an improvement over simpler 
approximations where they have the same value. However, the values of Nl and are 
rather different from the simulations, and this comes from ignoring higher correlations 
of the numbers of proteins and the state of the switch. Only when u^b ^ are the 
values in close agreement with the simulation values as expected in this limit where the 
mean field theory is exact. 

Mean field theory can, in principle, be improved by considering higher order 
moments and correlations. However, the algebra soon gets quite complicated. 

4. Exact perturbative solution 

In this section we develop a perturbative approach that allows the steady state 
probabilities Ps{Ni^ N2) to be computed as a power series in the unbinding rate u. 

4.1. Formal solution 

We begin by considering the formal solution of the master equation system (1-3). To 
transform the system of equations into a system of partial differential equations, we take 
the generating function of the different probability distributions: 

00 00 

Ksizi, ^2) = E E Zi'z2' PsiNi.N,) (44) 
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where = 0, 1, 2. In this way we obtain a system of hnear partial differential equations 
with non-constant coefficents: 

g{z, + z,- 2)Ko + [d-id + b)z^] ^ + [d-{d + 6)^2] ^ 

OZi (JZ2 

+UZ1K1 + UZ2K2 = (45) 
[g^,^ _l)_u]K, + d{l - zi)^ + d{l - Z2)^ + = (46) 

OZi OZ2 OZ\ 

[g{z2-l)-u]K2 + d{l-z^)^^d{l-z^)^ + h^ =0 (47) 

OZi UZ2 UZ2 

where the right-hand side terms have been set to 0, since the stationary probabilities 
Ps{Ni^ N2) are the main quantities to be determined in this paper. 

The second and the third equation of the system work in a completely analogous 
way, because of the symmetry of species 1 and 2, so it will be enough to deal with 
the first two equations. We also note the symmetries ^^2(^1,^2) = ^1(^2,^1) and 

^^0(^1,^2) = Kq{z2,Zi). 

In appendix A we give a formal solution to the system (45-47). However, in practice 
it is not clear how to actually compute e.g. probability distributions from this solution. 
In order to do this we develop instead a perturbative approach. 

4-2. Perturbative approach 

In this section we develop a perturbative approach to the problem of finding the exact 
stationary state. To do so we require a suitable small parameter of the model which we 
choose to be the unbinding parameter. In the 2/^0, the exact solution is simple: if, 
for example, one protein of type 1 is bound, the proteins of this kind will obey the usual 
Poisson distribution regulated by death and birth terms, while the number of proteins 
of type 2 will just decay to 0. This limit is the starting point for a perturbative solution, 

wherein the probability distribution will be expanded in a power series of w. 

00 

^5 = E""^"^- (48) 

n=0 

Owing to the symmetry of the system ^2(^1, ^2) = A (^2, ^1) we need only consider 
Po and Pi. 

Writing out the expansion explicitly we have 
Pi = + nP^'^ ■ ■ ■ 

Po = uP^'Hu^P^'\.. ^^^^ 

Note that the constant term P^^ = since Pq = in the limit of no unbinding. 

This approach also makes sense when the typical E. coli values for the parameters 
are considered [7]: 

^ = 0.05 rf = 0.005 6 = 0.1 2i = 0.005 (50) 

u is, along with rf, the smallest of the parameters. An expansion in 1/6 could also be 
developed, but we find the expansion in u more convenient. 
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4-3. Zeroth order 

In the zeroth order of the u expansion of the stationary master equation (3) (with l.h.s 
set to zero) we find 

= ^[P°(iVi - 1, N2) - P^{Ni,N2)] + d[{Ni + l)P°(iVi + 1, N2) + {N2 + l)Pr(iVi, N2 + 1) 
-{Ni + N2)P^{N,,N2)] . (51) 
We define the generating function 

CO CO 

k['\z„ .2) = E E ^^'^'PS'HNr, N,) (52) 

Ari=0 iV2=0 

which obeys 

= giz^ - 1)^^ + dil - zi)^ + d{l - Z2)^ (53) 

the solution of which is independent of Z2'. 

Kf^ =ciexp^zi (54) 

where Ci is a constant to be determined. Analogously K^^ {z) = 02 exp |z2- Applying the 
normalization condition K^\l) + K^\l) = 1 and the symmetry consideration Ci = C2 
leads to 

k[^^ = - exp exp ^zi . (55) 
Expanding as a power series in Zi , Z2 yields 



PiiNi,N2) = -^;^[^) 5iv,o. (56) 

This is a Poisson distribution for A^i with mean g/d^ with N2 fixed to be zero. The 
normalisation factor 1/2 is so that P(°)(iVi, A^s) = Pi^\Ni,N2) + P^^\Ni,N2) is 
normalised to unity. 

4^4' General Formulation 

We substitute the expansion (48) into the stationary master system (1-3) with time 
derivatives set equal to zero. Arranging orders of u the equations may be written as 

CsP^^\n„N2) = -f^''\N„N2) (57) 

for S = 0^1 where the action of the linear operators Cs is 

CoPt\N,,N2) = 9[Pt\N, - 1, N2) + Pt\N,, iV2 - 1) - 2Pt\N,, N^)] 



+ d 



(iVi + l)Pt\N, + l,iV2) + {N2 + l)Pt\Ni,N2 + 1) 



{N, + N2)Pt\Ni,N2) -b{Ni + N2)Pr{Ni,N2) (58) 



in). 



jC,Pt\N,, N2) = g[Pt\Ni - 1, N2) - Pt\N„ N2)] 



+ d 



{Ni + l)Pi''\Ni + 1, A^2) + {N2 + l)Pt\Ni, N2 + 1) 
>(»)/ 



-iN^ + N2)Pr'iN^,N2) (59) 
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and the inhomogenous terms are 

flr\N^, N2) = Pf''-'\N, - 1, N,) + Pt-'\N^, N2 - 1) (60) 
f[-\N,, N,) = - Pi""-'^ + b{Nr + l)Pt\N, + 1, N2) . (61) 

As noted above, we can first determine the zeroth order Pi^\Ni^ N2) and P2^\Ni^ N2) = 
pI^\N2^ Ni) as functions of the parameters of the model. Then, owing to the form of 
the equation (60), /q"^^ is determined. This allows us to solve for Pq"^^ which in turn 
determines f^^"^ and allows us to solve for P^^K Continuing in this fashion the rest of 
the probabilities will be found following the structure: 

) ^ Po^^^ ^ Pi'^ ^ Po^^) ^ Pf ) ■ • • (62) 

In general, Pq^^ will be found before P^^\ 

Having laid out the general perturbation scheme and established the zeroth order 
(n = 0) solution, we now outline how equations (57) can be solved. We note that for 
each switch state S the same linear operator Cs appears at all orders n. This means 
that the homogeneous parts of the equations are independent of the order (with only 
the inhomogeneous term on the right hand side varying between orders) and that they 
only have to be solved once. 

4-5. Green function for Cq 

Let us define a Green Function Qo{Ni^ N2\N^^ N2) for the operator Cq through 

CoQoiN,, N^lNl TVO) = -5^^,^o5^^,^o (63) 

SO that the solution of (57) may be written 

Pt\NuN2) = J2 fi''\N',^N',)Q,{N,^N2\N',^N',) . (64) 

We define a generating function 



00 00 



Ko(^i,Z2|iV°,iV°) =J2Y1 ^^^2^^Qo(iVi,iV2|iV°,iV°) (65) 

Ni=0N2=0 

the equation for which is obtained by summing (63) 

a{zi)^^ + a{z2)^^ + g{zi + Z2 - 2)Ko = -z^' z^^ (66) 

UZi OZ2 

with a{zi) = d — {d + h)zi. 

In order to solve (66) we use the method of characteristics (see e.g. [27]). The 
characteristic equations are 

^^«(..) (67) 

^ = «(Z2) (68) 

as 

^ = -,(., + . 2 -2)i^o-.f°.f (69) 
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where 5 is a time-like parameter and (67,68) define characteristic curves along which 
the partial differential equation (66) reduces to the ordinary differential equation (69). 
Solving (67,68) yields the curves 

where 

V = e-(^+^)^ (71) 
and A^B are two constants to be fixed. 



Equation (69) can be rewritten as an ordinary differential equation in v as 



d 



2pb g(A+B)^ 



1 2^^,-1 _9iA±B)_ f d . \ ^ ( d 



^ ^ d + b \d + b J \d + b ' ^ ^ 



^1 / rl \ ^2 



d^ 

To integrate (72) we choose an end-point of the integration as = 1 and set this to 
correspond to an arbitrary point in the Z1-Z2 plane. This fixes the two constants A^B 
as 

d ^ d , , 



Thus we obtain 



1 2^7b -. ( q f 2d 

Mz,, = I dvv^- exp (1 - .) 



{l-v) + vzi] (-—-{l-v) + vz2] . (74) 



d + b^ ' ) \d^h 
We now expand as a power series in z\, Z2 to obtain Qo{Ni, A^2|-^i , ) from (65) 

«,(^..;V,<^»)^l^X^i(j^) {^) (75) 



p=0 

T- ' 



rl \d+bj \N9 - rj \ d + b 

r=0 



where /(a, /3; x) is defined as the integral 

Iia,l5;x)= f dvv"-\l-vf-'e^^ = ^^P^ ,F,ia,a + /3;x) , (76) 
Jo 1 (a + p) 

and r(2;) and iFi{a,b;z) are the usual Euler Gamma function and confluent 
hypergeometric function respectively. 
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4.6. Green function for C\ 

If we define a Green function for the operator C\ through 

>Cigi(iVi, Nl) = -5^,,,^o5^^,,^^o , (77) 

this equation only has a solution when A^2 7^ 0. To see this we note that summing the 
left hand side of (77) over all A^i,A^2 yields zero i.e. the operator conserves probability. 
Therefore one cannot solve (77) or indeed (57) for an arbitrary right hand side; one 
requires that the sum of the right hand side of (57) over all A^i,A^2 yields zero. However, 
since the null space of the operator Ci is concentrated on A^2 = (i.e. the stationary 
state in (56) is proportional to S^^^q) we can find a solution of (77) for N2 > 0. 

It is simplest to proceed by considering the solution Pi for an arbitrary right hand 
side -h{Ni,N2) 

CiPi{N,,N2) = -h{N,,N2) (78) 

that satisfies 

00 00 

Y,J2HNi,N2) = 0. (79) 

Ni=0N2=0 

We define generating functions 

00 00 

Kr{zr,Z2)= H 'i' Pi (Ni , N,) (80) 

Ni=0N2=0 
00 00 

H{z,,Z2) = J2 J2''i'''2'h{N,,N,) (81) 

Ni=0N2=0 

then summing (78) yields 

d{l -z,)^ + d{l -Z2)^ + g{z, - l)Ko = -H{z,, z^) . (82) 

In order to solve (82) we again use the method of characteristics. The characteristic 
equations are this time 

^ = rf(l - ^i) (83) 
-T^ = c^(l - ^2) (84) 

A 

= -g{zi - 1)K, - H{zu Z2) . (85) 

as 

Solving the first two equations yields 

zi = l + Av Z2 = l + Bv (86) 

where now 

V = e-^^ (87) 

and B are constants to be fixed. The final equation becomes 
d 



d^ 



H{zi, Z2) 
dv 
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We choose the integration to be from ^ = to = 1 where v = corresponds to 
Zi = Z2 = I and v = 1 corresponds to an arbitrary point in the Z1-Z2 plane which 
imphes A — z\ — \ and 5 = ^2 — 1- We then obtain the solution of (88) 

^a \ (^l-l)f(l-^) 

a Jq V 

Expanding as a power series in ^1,^2 implies 

P.iNuN^) = K{1, l)e-^Mi^M^5r,,,o 

d \N2j^\dJ r\\Ni-r 

X I{Ni + N2 - r,p + q - Ni - N2 + 2r + 1; g/d) . (90) 

As discussed above, for A^2 > we can define the Green function (77) by means of which 
the solution of (78) may be written 

A(iVi, N2) = 5] gi(iVi, N^Wl N^)h{Nl N^) . (91) 
Then we can read off the Green function from (90) as 

X I{Ni + N2-r, + - Ni - N2 + 2r + 1; g/d) 

For A^2 = the solution (90) reads 

Pi(iVi,0) = K(l,l)e-«/^- 



,ig/dr^ 



-g/d 



2^2^ d 2^\d) r\\Ni-rJ 



p=0 q=0 r=0 

X I{Ni-r,p + q- Ni + 2r + l;g/d) . (93) 

Some care is required with the r = Ni term of the sum in (93) since the integral 
/(0,/3;x) does not converge. However the property (79) of the function h implies that 
the coefficient of the offending integral is zero. To see this one can write the r = Ni 
term of (93) as 



oo oo 



i:f^(f)"'^/-e-5-.„-.,:-.r-. 

p=0 q=0 ^' s=l ^ / 
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In the final equality, the binomial expansion of (1 — vY^^ has been used with the term 
5 = not present since its coefficient vanishes due to (79). All the integrals in (94) then 
converge. 



4-7. First- order results 

The first-order contribution to the stationary probability Pq^^(A^i, A^2) is given by: 

Pi'\N,,N2)= fi'\NlN'M{N,,N,\N',,N',) (95) 



with 



/r«iV2") = ^exp(-^ 



(iVO - 1)! 



iV"-l 



(iVO - 1)! 



.(96) 



Although Qo{Ni^ A^2|A^i , ^^2) expressed in terms of known integrals in (75), it is 
advisable to go back to the explicit expression of the integrals (76) to evaluate the sums 
appearing in (95) more easily. In that way, Pq'^^(A^i, A^2) can be written as: 



_9__J2gd_ 



le (d+b)2 ^ 1 



2 d+b 



N^=0 
N2 



d 



Ni 



V- 



9 



d + b J \Ni - pi \d + b 



d 



X 



2gb 



1 _ yf?-Ni+N,+2p^^v ^ ^^^^ 



where the label symm refers to the fact that there will be another term equal to the 
written one, apart from a switch in the variables and A^2 • 

In Appendix B it is shown how the expression may be simplified to the result 

N2 Ni . X Ni-m 



P^'\NuN2)= -exp 



1 



gjb - d) _ g_ 
[d + bf d 



1 



9 



dJ d + bN2\\d + b 



m=0 



9 



X 



1 



2gb 



d+b 
(d-b) 



d. 



{Ni - m)lml [d + b \{d + b)^ {d + bf 



+ ml 



'^gb , 9{d-b) 

^ +m,N2 + Ni-m + 1-^ ' 



symm . 



Xd + by ' ' ' ' {d + by 

Once we have Po^\ we can plug this result into the P/^'' equation which becomes 
for A^2 > 

Pi'\N,,N2)= A'\N'i,N^)Qi{Ni,N2\NlN',) (99) 

N°,N° 

where ^(iVi, N2) = -PI'\Ni, N2) + 6(iVi + l)Pi'\N, + 1, N2). 
We consider separately the two terms of /^[^^ . The first is 



(98) 



(100) 



2 \dJ 

Since the calculation with the Green function (see section 4.6) is only for A^2 7^ 
this term does not contribute and the only contribution comes from the second 
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term in f[^^ involving Pq^\ The resulting expression, obtained by substituting h 
b{Ni + 1)P^^\Ni + 1, A^2) in (90) with P^^^ given by (98), is 

\ OO OO ^ 

ATA - JL I yS"^^^ _ ^\ ^(7V0 + l)j 

1 1 / ^ \ '^2 / n \ ^1+1-^ 



(d + b) \d + b 



9 



2g_ 
d 

iVU=OiV^=0 

E 



(101) 



m=0 



+ ml 



I 



2gb 



(d + b)^ 
2gb 



m, + -m + 3; 



{d + by 



+ m, + A^i - m + 2; 



^(^-^) ' 
(rf + 6)2 
(d-b) 



1 — m)!m! 



(d + 6)2 



s?/mm(A^° + l,A^°) 



ox Wi 



H U i JW + iV2-r,iV° + iV°-iV,-iV, + 2r + l;^/ci) 



Where the symmetric term in this case corresponds to the term inside the curly brackets, 
exchanging the places of + 1 and coming from the symmetric term in Pq . For 
A^2 = we obtain the result shown in (93) and cannot be simplified further. 

Equations (93), (98) and (101) are the main results of this section and give closed 
form expressions for the first-order contributions to the stationary proabilities. The 
normalization constant K{1^ 1) appearing in (93) is obtained from the condition: 



. 



(102) 



E E [Pi'\NuN,) + Pi'\N,,N,) + Pi'\N,,N,) 

Ni=0N2=0 

Once the values of the probabilities Pj^^\Ni^ N2) have been obtained numerically, 
they are multiplied by the unbinding parameter u and added to the zeroth order 
probabilities. In that way, the probability distributions P^(A^i,A^2) can be computed 
and plotted up to the first order of the expansion. Figure 4 shows the probability 
distributions for different values. The agreement is visually good in the first two 
examples: typical E. coli values (4) and another interesting case, with smaller g. 
The order of magnitude is well reproduced in almost every point of the probability 
distribution. However, the approximation does not work accurately if the value of the 
unbinding parameter u is of the order or bigger than the other parameters. 

Although the first two examples of figure 4 seem visually in good agreement with 
the simulations, the best way to check this is to plot different slices of the probability 
distribution, that is the probability distribution of N2 where A^i is held constant (figure 
5). We choose the values of A^i to correspond to the slices with large probability mass, 
i.e. A^i = 0, 1,2 for E. coli values. Figure 5) shows that along the slice with greatest 
probability mass (A^i = 0) the analytical and simulation plots show good agreement. 
For A^i > 0, there is reasonable agreement for the E. coli values (figure 5 b)), whereas 
for larger u the agreement is not so good (figure 5 d)). 

Since the proposed method is general, calculations can be performed up to any 
necessary order to get better results. For example, in the case of E. coli in axes A^i = 1, 
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A^2 = 1 is enough to compute the second order to have very good results (figure 6). In 
general, the method can be iterated as many times as required. 



Analytical E. CoH Simulation E.Coli 




a) E. coli values. 

Analytical g=0.0 1 r.,u=0.000ri Simidalion g=0.015,u=0.0005 




h) g = 0.015, u = 0.0005, d = 0.005, b = 0.1. 

Analytical u =0.0r. Simulation u=0.05 




c) E. coli values with u = 0.05. 

Figure 4. Comparison between the analytical first order and the simulation 
distributions for different values of the parameters a) E. coli values (4). The probability 
distributions have the same shape and the order of magnitude is well reproduced in 
almost every point, b) ^ = 0.015, u = 0.0005, d = 0.005, b = 0.1 The order of magnitude 
is again well reproduced. In this case, the two peaks get closer to the origin as the 
ratio g/d is smaller, c) E. coli values with u = 0.05. The approximation at first order 
is no longer accurate as the value of u is no longer small compared to the rest of the 
parameters of the model. 

5. Conclusion 

The exclusive genetic switch presented in this paper represents a minimal model of two 
populations that compete in an indirect way, that is, through a common promoter site 
that when occupied by a member of one population stops the production of the other. 
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X X Analytical 
O O Simulation 



O o Analytical PCl.N^) 

X X Simulation P(1,N2) 

Analytical P(2,N2) 
+ ....+ Simulation P(2,N,) 



a) E. coli values. Slice A^i = 0. 



O O Analytical 
X X Simulation 



c) ^ = 0.015, u = 0.0005, d = 0.005, b 
0.1. Slice A^i = 0. 



b) E. coli values. Slices A^i = 1,2. 



O O Analytical 
X X Simulation 



d) E. coli values with u = 0.05. Slice 
N2 = l. 



Figure 5. Comparison between the analytical probability distributions P{N^ ^ N2) 
where is fixed and chosen to correspond to slices with largest probability mass. 
The agreement in slice Ni = is good for both figures a) and c). There is some 
quantitative difference in b), which corresponds to the E. coli values in slices with less 
probability mass. This is improved with second order calculations, as can be seen in 
figure 6. In d) the difference is clear, since the first order approximation is no longer 
accurate, as discussed in figure 4. The error in the values is negligible and, in all the 
cases, is smaller than the size of the used symbols. 



It is a non-equilibrium system because the microscopic processes are irreversible and 
detailed balance does not hold. The non equilibrium stationary state exhibits interesting 
properties of bistability i.e. generically at any given time the system is dominated by 
one of the populations. We have developed two main analytical approaches to study the 
stationary state: a mean field theory and an exact perturbative approach. In addition 
we have presented some exactly solvable limits. We also have studied the system and 
checked the analytical results by using Monte Carlo simulations. 

These simulations show that the symmetry of the system is always broken, that 
is, the system is always functioning as a switch, with two opposite states in which 
one population is much more abundant than the other. This holds generically except 
in the limit u^b ^ oo where the state is symmetric. As the relevant parameters 
for the switch become smaller, the probability of finding bound proteins decreases, 
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Figure 6. Comparison between the probability slices P(1,7V2) and P(2,A^2) from 
second order analytical calculations and simulations, for E. coli values. The second 
order is clearly enough to get accurate results. 

0.003 I ^ \ ^ \ ^ — 
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0.0015 
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but the distributions always remain asymmetric, showing no phase transition between 
asymmetric and symmetric states 

The mean field theory of section 3 is based on exact moment equations which are 
factorised at the level of second moments. The approximation scheme is constructed 
to be exact in the limiting parameter case of vanishing binding/unbinding rates 6, u 
with the binding constant k = ^ held constant The theory shows good agreement 
with some quantities obtained from the simulation, but there are some other quantities 
whose agreement varies and depends on the value of the parameters. The mean- field 
theory could be improved upon by systematically considering higher order moments and 
correlations. 

As a first attempt, to our knowledge, at an exact solution of this nonequilibrium 
system we have developed an exact perturbative approach which consists of an expansion 
in the unbinding rate u. We computed the expansion to first order and this has allowed 
us to obtain the whole probability distribution for the typical E. coli and other values, 
in good agreement with the simulations. Higher orders can be obtained systematically 
by iterating the proposed method. The analytical expressions at first order (98,101) 
already illustrate the complexity of the nonequilibrium stationary state. In principle, 
the Green functions that have been calculated in Section 4 allow the expansion to be 
carried out to arbitrary order although analytical expressions for the higher order terms 
in the expansions might be long or difficult to simplify. However, they still can be 
computed numerically by programming the proposed operations and, in principle, the 
method can be iterated as many times as necessary to get more accurate results. 

In this paper we have not attempted to study the dynamics, but it would be of 
interest to do so. In particular the dependence of the flip time or first passage time (the 
time for the system to change from being dominated by one population to the other) on 
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the parameter values is of interest. It may be possible to extend our analytical solution 
to consider such flip times and to provide estimate of flrst passage times between the 
two bistable states. It might be that the Green functions we have computed in section 
4 contain some useful dynamical information. 

The techniques that we have developed should be applicable to the understanding of 
other related systems and properties of genetic switches in general. The factorisation of 
the moment equation hierarchy is straightforward to implement to obtain a mean fleld 
theory; the exact perturbative approach is more involved but can be used in general 
for problems with similar probability distributions, as long as the Green functions are 
analytically solvable. Thus the techniques represent standard procedures to analyze this 
class of systems. 
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Appendix A. Analytical method for systems of linear PDEs. 



According to [27], there is a way to solve, or at least simplify, a system of linear first- 
order partial differential equations. First of all, considering a three dimensional space 
with a vector Z2) for the three probabilites, the system has to be written as: 

where A,B and C are 3x3 matrices, and d and K are column matrices. The requirement 
to solve the system is that detA{zi, Z2) 7^ or detB{zi, Z2) 7^ 0. 
In this case the expression of the matrices is: 

/ d-{d + h)zi \ 

A{zi,Z2)= b d{l-zi) 





B{zi,Z2) 



C{zi,Z2) 



\ 

( d-{d + b)z2 

d{l - Z2) 



d{l-zi) j 








\ 



( -g{zi + Z2 


V 



d{l-Z2) J 

2) —uzi —UZ2 

-[g{zi-l)-u] 

-[9{Z2-1) 



U 



J 



d{zi,Z2)= (A.2) 

This is not the most general system that could be written with this notation, since the 
matrix d is zero, and A and B only depend on Zi, ^2, respectively. 

The requirement to solve the system is that detA{zi^ Z2) 7^ or detB{zi^ Z2) 7^ 0, 
which is fulfilled almost in every point. Multiplying by the matrix the matrix A 
is transformed into: 



A' = B-^A 



d-{d+b)zi 
d-{d+b)z2 
b 

d{l-Z2) 
-b[d-(d+b)zi 



d{l-Z2){d-{d+b)z2) 





(l-^l) 

1-Z2 








(l-^l) 

1-22 



(A.3) 
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whose eigenvalues, written as columns in a matrix R are: 

/ \ 



R 



d(zi—l)+bzi 
d-{d+b)z2 
d-{d+b)zi 

1 



(A.4) 



1 

VI 

This hyperbolic system can be transformed with elementary matrix operations: 

C = B-^C K = R{zi, Z2)v (A.5) 

The method then states that the components of v obey the following system of equations: 



CijVj 



(A.6) 



where Cij are the components of the matrix C once we have performed the 



transformation with R: C = R~^C'R. 

As can be seen from the last equation, we have uncoupled the derivative terms 
and, even if this system cannot be solved analytically, it is easier to deal with it 
computationally. The problem itself can be written as: 
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dvi 
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d{zi - Z2){d - (6 + d)z2) 



(A.7) 



d — {d + h)z2 

Although this method in principle solves exactly the system of partial differential equa- 
tions it appears a formidable task to actually integrate equations (A.7). 
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Appendix B. Derivation of expression (98) 

In this appendix we give the detailed derivation of (98). We begin from (97) 

y.1- {^-^^ dw^^^^-'-^l - vf°-"'+'-'+^'e^" + symm (B.l) 

where the label symm refers to the fact that there will be another term equal to the 
written one, apart from a switch in the variables and • 
Now defining: 

gb gd X d 1 — V Q 

civ)=v^^e^^^ \ m = Ni-p, r(v) = - — , s(v) = -^(l-v)(B.2) 

we arrive at the expression (for convenience we will drop the dependence of the previous 
functions on v): 



(B.3) 



^\mj {N^-m)\ {N^y. ^ 



m=0 

Separating the parts of the expression that can be summed, the following 
simplification can be obtained by changing the order of the sums appropriately: 



2^ [dJ (m9 - 1^1^ ' 2^ 



N°=0 



dJ (iV?-l)! ^,\mj {N,-m)\ 



Ni 



X -1 ^ y> + (B.4) 



rf/ (A^i — m)!m! n! 

m=0 ^ ^ n=0 



dJ - mV.ml ^ ^ 



where uj = vg/d. Note that this sum could be simplified further, but that the 
simplification will not allow us to perform the integration over v in closed form. In 
this and following equation, we will try to obtain the simplest expressions globally, 
knowing that simplifying one part can lead to further complications in another. 

Plugging this sum into the Pq^^ equation and writing all the explicit forms of the 
functions, we obtain the result (98). 



